# This code replicates Table 8
rm(list=ls())
setwd("")
library(foreign)
data <- read.dta("data.dta")
data$traveled <- data$traveled-1

# Define variables
treats <- c("treatment_catholic" , "treatment_adventist" , "treatment_mixed" , "treatment_peruana" , "treatment_maranatha")
vars <- c("age", "male", "education", "income", "income_none", "married", "relationship", "single", "farmer", "housewife", "internet", "traveled", "exp1_leave", "exp1_leave_state", "exp1_leave_church", "accepts_wrongdoings_leader", "exp4_swayed", "exp4_swayed_state", "exp4_swayed_church", "state_projects_count", "community_projects_count")

# Write helper function
get_stats <- function(var, treats, data) {
  c(mean(data[, var][data[, treats] ==1 ], na.rm = T), sd(data[, var][data[, treats] ==1 ], na.rm = T))
}

# Loop through variables
btable <- lapply(treats, function(y) t(sapply(vars, function(x) {
  get_stats(x, y, data = data)
})))

# Convert to data frame
btable <- data.frame(do.call(cbind, btable))

# Adapt percentages
perc <- c("male", "income_none", "married", "relationship", "single", "farmer", "housewife", "internet", "traveled","exp4_swayed", "exp4_swayed_state", "exp4_swayed_church", "exp1_leave", "exp1_leave_state", "exp1_leave_church", "accepts_wrongdoings_leader")
btable[rownames(btable) %in% perc, ] <- btable[rownames(btable) %in% perc, ] * 100

# Add varnames
btable$variable <- rownames(btable)

# Add N
n <- sapply(vars, function(x) sum(!is.na(data[, x])))
btable <- cbind(n, btable)

# Change colnames
cn <- paste0(rep(treats, each = 2), rep(c("_mean", "_sd"), 2))
colnames(btable)[c(-1, -ncol(btable))] <- cn
btable <- data.frame(round(btable[, -which(colnames(btable) == "variable")], 1))
rownames(btable) <- gsub("_", " ", rownames(btable))

# Output
library(xtable)
xtable(btable,digits=1)
